#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/* http://gmplib.org/*/
#include <gmp.h> 

/*IMPORTANT NOTE*/
/* Compile with the line: gcc calculo_pi.c -o calculo_pi -lm -lgmp*/	

static const int terms = 50;
int x = 1;

//modify this value
long int work_unit = 3000;

long int units_to_complete;


/* prototypes */
void get_pi();
void factorial(unsigned long int n, mpz_t fact);


void factorial(unsigned long int n, mpz_t fact)
{
	unsigned long int  num = n;
	mpz_fac_ui (fact, num);
}

void get_pi()
{	

	units_to_complete = terms*work_unit;
	long double pi_sum = 0.0;
	long double partial_pi;
    long int n;
    
	/*variables required for the operation*/
	unsigned long int two;
	unsigned long int n2;
	unsigned long int n2p;
	two = 2;
	
    /* Number of terms to complete! */    
    for(n = 0; n<units_to_complete; n++){
        
		
		/*setting up the objects required*/
		mpf_t dividend_float, divisor_float;
		mpf_init(dividend_float);
		mpf_init(divisor_float);

		mpz_t dividend;
		mpz_init (dividend);
	
		mpz_t base;
		mpz_init (base);
		mpz_set_ui(base, two);
	
		mpz_t divisor_A;
		mpz_init (divisor_A);
	
		mpz_t divisor_B;
		mpz_init (divisor_B);

		mpz_t divisor_C;
		mpz_init (divisor_C);
	
		mpz_t divisor_T;
		mpz_init (divisor_T);
	
		mpz_t divisor_F;
		mpz_init (divisor_F);
	
		mpf_t partial_result;
		mpf_init (partial_result);

		/* operations */
		n2 = 2*n;
		n2p = 2*n+1;
		
		factorial((2*n), dividend); 
		mpz_pow_ui(divisor_A, base, n2);
		factorial(n, divisor_B);
		mpz_pow_ui(divisor_B, divisor_B, two);
		mpz_set_ui(divisor_C, n2p);
		mpz_mul(divisor_T, divisor_A, divisor_B);
		mpz_mul(divisor_F, divisor_T, divisor_C);		
		mpf_set_z(dividend_float, dividend);
		mpf_set_z(divisor_float, divisor_F);
		mpf_div(partial_result, dividend_float, divisor_float);
		
		partial_pi = (long double)mpf_get_d(partial_result);
		
		pi_sum = pi_sum + (2*partial_pi);
		
		printf("-> Current Value of Pi: %Lf\n", pi_sum);
		
		
		/*clean space used by the objects*/
		mpz_clear(dividend);
		mpz_clear(base);
		mpz_clear(divisor_A);
		mpz_clear(divisor_B);
		mpz_clear(divisor_C);
		mpz_clear(divisor_T);
		mpz_clear(divisor_F);
		mpf_clear(partial_result);
		mpf_clear(dividend_float);
		mpf_clear(divisor_float);
	}
	printf("\n***** FINAL Value of Pi: %Lf *****\n", pi_sum);
	
} 

int main()
{
  
  printf("*** Operation to determine the value of Pi ***\n");
  printf("-> Number of work units: %li\n", work_unit);
  printf("-> Number of iterations: %lin", work_unit*terms);
  get_pi();
  return 0;
}

